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ABSTRACT 

Algorithms  for  solving  the  isotonic  regression  problem  in 
two  dimensions  are  difficult  to  implement  because  of  the 
large  number  of  lower  sets  present.  Here  a  new  algorithm  for 
solving  this  problem  based  on  a  simple 
proposed  and  shown  to  converge  to  the 
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1.  INTRODUCTION.  Algorithms  for  calculating  the  least 
squares  isotonic  regression  function  have  received  a  great  deal 
of  attention  in  the  literature  and  six  such  algorithms  are  dis¬ 
cussed  in  Section  2.3  of  Barlow,  Bartholomew,  Bremner  and 
Brunk  (1972).  In  situations  where  there  is  one  independent 
variable  all  of  the  algorithms  work  very  efficiently.  Perhaps 
the  most  widely  used  algorithm  is  the  "pool  adjacent  violators 
algorithm"  which  is  applicable  only  in  the  case  of  a  simple 
linear  ordering  or  an  amalgamation  of  simple  orderings.  In 
many  isotonic  regression  problems  we  have  more  than  one  inde¬ 
pendent  variable  present  and  are  concerned  with  partial  order¬ 
ings.  An  Important  example  involves  the  prediction  of  success 
in  college.  Usually,  this  prediction  is  based  upon  several 
independent  variables  such  as  rank  in  high  school  graduating 
class  and  score  on  a  standardized  examination  such  as  the  ACT 
composite  and  is  measured  In  terms  of  a  predicted  grade  point 
average  or  predicted  probability  of  obtaining  a  particular  GPA 
or  better.  The  predicted  value  Is  usually  obtained  by  regres¬ 
sion  methods  and  is  assumed  to  be  nondecreasing  In  each  inde¬ 
pendent  variable.  The  isotonic  regression  function  has  been 
found  to  compare  very  favorably  with  other  techniques  with 
respect  to  predictive  accuracy  (cf.  Perrin  and  Whitney  (1976) 
and  Kolen  and  Whitney  (1978)). 

Some  of  the  algorithms  described  in  Barlow  et  al.  are 
applicable  to  the  case  of  computing  the  doubly  nondecreasing 
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least  squares  regression  function  but  the  number  of  computa¬ 
tions  required  can  become  prohibitive.  For  example,  consider 
the  minimum  lower  sets  algorithm  described  in  Section  2.3  of 
Barlow  et  al.  Suppose  one  of  our  two  independent  variables  has 
a  possible  values  and  the  other  has  b  possible  values.  By 
counting  paths  from  the  upper  left  hand  corner  to  the  lower 
right  hand  corner  of  our  a  xb  grid,  it  follows  that  the  number 

of  lower  sets  is  equal  to  •  If  a  =  b  this  number  is 

—  1/2  a 

approximately  (an)  -H  by  Stirling's  formula.  Thus  if 

a  =  b  =  20,  and  if  consideration  of  each  lower  set  were  to 
require  one  microsecond  of  computer  time,  then  finding  the 
first  level  set  would  require  2312  minutes  or  38.5  hours  of  CPU 
time.  (One  microsecond  seems  conservative  in  light  of  the  fact 
that  computation  of  the  average  value  over  that  set  would  take 
at  least  two  multiplications,  two  additions  and  a  division  and 
the  comparison  would  require  a  subtraction.  The  present  stan¬ 
dard  for  making  such  predictions  is  four  arithmetic  operations 
per  microsecond.)  Moreover,  if  the  first  level  set  is  small 
(as  it  would  be  with  good  data)  the  second  cycle  is  nearly  as 
difficult  as  the  first,  etc. 

Since  the  doubly  nondecreasing  regression  function  is  so 
difficult  to  compute,  researchers  have  proposed  using  ad  hoc 
estimators  based  upon  one  dimensional  smoothings.  (The  number 
of  computations  required  for  one  dimensional  smoothings  is 
essentially  linear  in  the  number  of  entries.)  Makowski  (197*0 
studied  consistency  properties  of  estimators  obtained  by 

p 

r 
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successive  one  dimensional  smoothings.  Kolen,  Smith  and  Whitney 
(1977),  Perrin  and  Whitney  (1976),  and  Kolen  and  Whitney  (1978) 
proposed  two  different  techniques  for  producing  estimates  which 
are  nondecreasing  in  each  variable.  One  of  their  techniques 
was  to  first  do  one  dimensional  row  smoothings.  After  all  rows 
had  been  adjusted,  reversals  in  the  columns  were  adjusted  by 
the  same  method.  They  then  returned  to  the  original  table  and 
did  one  dimensional  column  smoothing  followed  by  row  smoothings. 
Neither  smoothing  necessarily  produces  a  doubly  nondecreasing 
table  so  they  averaged  the  two  results.  (The  average  is 

not  necessarily  doubly  nondecreasing  but  was  for  their  data.) 
This  method  was  applied  to  the  problem  of  estimating  the  prob¬ 
ability  of  obtaining  a  "B  or  better"  GPA  for  entering  college 
students.  The  data  is  presented  in  Table  1.  The  two  entries 
are  the  total  cell  frequencies  and  the  observed  relative  fre¬ 
quencies.  We  note  that  there  are  a  number  of  "reversals," 
even  with  a  relatively  large  sample  size.  The  smoothed  esti¬ 
mates,  by  their  method,  are  presented  in  Table  2  and  the  iso¬ 
tonic  regression  function  with  weights  equal  to  frequencies 
in  Table  3-  Note  that  not  only  the  estimates  but  also  the  level 
sets  are  different.  These  level  sets  are  very  useful  for  making 
inferences  about  equivalent  scores  within  the  table. 

Several  of  the  algorithms  discussed  in  Section  2.3  of 
Barlow  et  al.  are  basically  methods  of  finding  linear  orders 
which  are  consistent  with  a  partial  order.  We  have  not  been 
able  to  write  a  program  which  implements  any  of  these  methods 


TABLE  1 


The  probability  of  making  a  "B  or  better"  GPA 
(top  number  =  total  cell  frequency;  bottom  number  =  relative  frequency) 


ACT 

High  School 

Grade 

Point  Average 

Composite 

0  to  1.55 

1.56  to  2.25  2.26  to  2 

1.95  2.96  to  3.65 

3.66  to  U.00 

+ 

CO 

CM 

0 

7 

10 

1*7 

1*1* 

.0000 

.2857 

.2000 

.571*5 

.8861* 

23-27 

7 

56 

88 

180 

81* 

.0000 

.1250 

.1818 

.2833 

.5238 

18-22 

23 

166 

152 

ll*9 

33 

•  OU35 

.0301 

.0721* 

.191*6 

.1212 

13-17 

27 

1U9 

96 

6l 

1* 

.0000 

.01*70 

.0313 

.01*92 

.5000 

0-12 

10 

.0000 

57 

.0000 

33 

.0606 

7 

.0000 

0 

.0000 

TABLE  2 

The  probability  of  making  a 
Kolen  and  Whitney 

"B  or  better 
Method 

"  GPA 

.0311* 

.2353 

.2353 

.571*5 

.8861* 

.0311* 

.1250 

.1818 

.2833 

.5238 

.0311* 

.0375 

.0721* 

.1867 

.1931* 

.0375 

.01*02 

.01+93 

.1781* 

.0383 

.01*21 

.01*25 

TABLE  3 

The  probability  of  making  a  "B  or  better"  GPA 
Least  squares  isotonic  regression  (weights  =  cell  frequencies) 


for  a  doubly  nondecreasing  regression  function  in  a  reasonable 
amount  of  time. 

In  this  paper  we  present  an  algorithm  for  calculating  the 
least  squares  isotonic  regression  function  which  is  increasing 
in  each  of  two  variables.  This  algorithm  uses  successive  one 
dimensional  smoothings  and  is  very  efficient  and  easy  to  program. 
To  illustrate,  we  applied  the  algorithm  to  a  20  by  20  table  of 
random  numbers.  The  isotonic  regression  was  obtained,  correct 
to  four  significant  digits,  after  400  iterations,  required  39 
seconds  of  CPU  time  at  a  cost  of  one  dollar  and  thirty-five 
cents.  This  algorithm  is  described  in  Section  3.  In  Section  2 
we  summarize  some  well  known  properties  of  isotonic  regression 
which  will  be  used  in  the  proof  that  the  algorithm  yields  the 
desired  result. 


2.  SOME  PRELIMINARIES.  We  let  0=  £ ( 1 , J ) ;  i=l,2,--*,aj 
j  *1,2,* ••,b]  and  define  the  partial  order  «  on  0  by 
(i,j)  «  (k,^)  if  and  only  if  i-k  £  0  and  j - ^  £  0.  We 
denote  an  arbitrary  real  function  whose  domain  is  0  as  a 
matrix ,  i . e . , 


G  =  (g±1 )  =  (g( (i » J ) ) ) ,  1  = 1,2, • • • ,a;  J  = 1 , 2 ,  • • •  ,b . 

(Note  that  this  is  not  the  usual  matrix  notation  where  g^ 
refers  to  the  entry  in  the  i—  row  and  j  — 


column . ) 


We  say  that  a  function  P  :  fl  — »  R  is  isotonic  or  order  preserv¬ 
ing  if  (i,j)  «  (k,i)  implies  fij  s  This  is  equiva¬ 

lent  to  requiring  that  P  be  nondecreasing  along  both  rows 
and  columns.  The  least  squares  isotonic  regression  problem  is 
to 

Minimize  I*  ) 2  "tJ 


for  P  belonging  to  the  class,  K,  of  Isotonic  functions 


where  w^  >  0  and  G  are  given. 

Since  the  class  of  isotonic  functions  forms  a  closed  convex 
cone,  it  is  well  known  (for  example,  see  Theorem  7.8  in  Barlow 
et  al.)  that  the  solution  to  the  isotonic  regression  problem. 


say  G*, 

is  characterized  by  the 

(2.1) 

(2.2) 

Ll,J(eU  -<!)h 

w, 


=  0, 


and 


w,  ,  SO, 


for  all  functions  H  €  K. 
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3.  THE  ALGORITHM.  The  algorithm  which  we  propose  requires 
only  the  ability  to  solve  the  isotonic  regression  problem  with 
the  usual  nondecreasing  order  (in  one  dimension)  along  rows  and 
columns.  Our  algorithm  is  given  as  follows: 


1.  Let  denote  the  isotonic  regression  solu- 

A(  i ) 

tion  of  G  =  (g^  . )  over  rows,  i.e.,  G  minimizes 


UU  subject  to 


s  f  ,  for 
aj 


j  =  !,•••, b.  We  call  R(1)  =  (r  (1))  =  (gij1'>  ~  gi  j  } 


the 


I 

i 


first  set  of  "row  increments." 
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Let  =  (g^j 1  )  denote  the  isotonic  regression  solu¬ 
tion  over  columns  from  the  initial  values  G+R^^,  i.e., 

minimizes  L  (g^j  tr^1^  “  ^ij  ^wij  subject  to 

f il <  fi2  £ • • •  £ fib  for  i  = 1 , • • • ,a.  We  call 

-(G+R^)  the  first  set  of  "column  increments." 
Note  that  G ^ 1  ^  =  G  +R(1^  +C^1). 


th  n ) 

3.  Etc.  At  the  beginning  of  the  n  cycle,  we  obtain  G 

by  isotonizing  G+C^n_1')  over  rows.  The  nth  set  of  row 
increments  is  defined  by  R^n^  =  G^n^  -(G+C^n  so  that 

=  G+C^n-1^  +R^n\  We  then  obtain  by  isotoniz¬ 

ing  G+R(n)  over  columns.  The  nth  set  of  column  incre¬ 
ments  is  given  by  C(n)  =  -(G+R^),  or  equivalently 

G^  =  G+R(‘n'>  +  C^n\ 


The  utility  of  the  algorithm  lies  in  the  following  theorem. 

Theorem  3.1.  Both  G^  and  G^  converge  to  the  true  solu- 

* 

tion  G  as  n  -*  00 . 


Proof.  If  we  denote  the  inner  product  norm  as 


l[F|| 


1 

(F,P)2 


a  b 

(  Z  Z 
i=l  J-l 


y 


we  first  show 
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(3-D  li£(n)||2  *  l!3(n)U2  *  ||&(n+1)  ||2  for  all  n. 

To  establish  some  additional  notation,  we  denote  the  "row  cone" 

by 

Kr  -  'F'  flJ  £f2j  5faj  f°r  3  ’I. •••.»], 

and  the  "column  cone"  by 


K 

c 


Cp;  rn*f12 


s  f  for  i  = 1 ,  •  •  ■  ,a) . 


The  respective  dual  cones,  as  discussed  in  Barlow  and  Brunk 
(1972),  are 


K*w 

r 


(H; 


1f1hijfijwij 


£  0 


for  j  =  1 ,  •  •  •  ,  b  ;  for  every  F  6  ] 


and 

K»w 

c 


b 

£ H ;  X  h .  .  f  w  .  £  0  for  i  =  1 ,  •  •  •  ,a ;  for  every  F  £  K  ). 


Since  -R^n+^  is  the  projection  of  G  +C^n^  onto  K*w,  the 
work  of  Barlow  and  Brunk  guarantees  that 


-R(n+1)  minimizes  ||  G  +  C(n)  -  F [| 2  for  F€K*W. 
Similarly , 

-C(n)  minimizes  |JG  +  R(n)  -  F |J 2  for  FCK*W. 
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Therefore,  since  -R(n)  6  K*w  and  -C(n_1)  €  K*w, 

*  r  c  * 

II Q  +R(n)  -  (-C(n-1))|j2  S  || G  +R(n)  -  (-c(n)  )  II2  i  II G  +  C(n)  -  (-R(n+1)  )  II2 

which  is  equivalent  to  (3.1). 

Next  we  show  that  { C ^ n ^ }  and  [R^n^]  are  bounded.  If 

not,  let  (ig,j  )  be  a  minimal  point  in  0  (with  respect  to 

our  partial  ordering  <<)  such  that  either  {r:  .  ]  or 

10*']0 

{cfn\  3  is  unbounded.  Say  there  exists  a  subsequence  {n.  ] 

VJo  *o 

such  that  r.  — ♦-«>.  (Since  E  r^n]  £  0  for  all 

10,J0  i=l  >J0 

(n.,  ) 

n  (see  Barlow  and  Brunk  (1972)),  r.  — ► 00  would  contradict 

10,J0 

the  fact  that  (i0,jQ)  is  minimal.)  But  this,  together  with 

G^n^  =  G+R^n^  +C^n^  and  the  fact  that  G^n^  is  bounded  in 

(n  ) 

norm  (cf.  (3-D)  implies  that  C.  .  — »  °°.  This,  In  turn, 

10,J0 

contradicts  the  fact  that  (1q,Jq)  is  minimal  since 
0 

E  cfn\  s  0  for  all  n. 
j=l  0’J 

Projections  onto  convex  sets  are  distance  reducing  so  that 
||C  (i  )  -  C  (i_1 )  ||2  =  || G  +C(1)  -  (G+C(1-1)  )  ||2 
(3.2)  *  ||R(1+1)  -R(1)||2=  || G  +R(i+1)  -  (G+R(i))||2 

s  ||C(i+1)  -C(1)||2  for  all  i. 


We  now  show  that 
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(3-3)  |iR(i+1)  -R(1^|l2  (and  hence  ||c(i+1)  -C(i)||2)-*  0 

as  i  — » 

If  (3-3)  were  not  the  case,  there  would  exist  ( i Q , j Q )  €  fi 
and  e  >  0  such  that 

(3.*0  |r.^i+:P  -r.^,  |  >  e  for  infinitely  many  i. 

X0’J0  10,J0 

However,  since  £ R ^ n ^ }  is  bounded,  there  exists  a  finite  M 
such  that 

(3.5)  |r f 1 }  .  -r.(J).  |  <  M  for  all  i,J. 

0,J0  0’J0 

If  we  write 

(3.6)  ||R(i+1)  -  R(1)  !|2  -  !|C(1+1)  -C(i)|j2 

=  ||G  +R(i)  +c(1)  -  (G+R(i+1)+C(1+1)  II2 

+  P(G+R^+C^-(G+R^1+1^  +  C^1+1^ )  ,C^i+1^  -C(1)), 

the  left  side  of  (3.6)  converges  to  0  since  by  (3-2)  both 
terms  converge  to  the  same  quantity.  The  last  term  of  the  right 
side  is  nonnegative  by  (2.1)  and  (2.2).  Thus 

(3-7)  (R^i+1^  -R^)  +  (C(i+1)  -C(1))  — *  0  as  i  -!-♦  •. 

In  similar  fashion,  beginning  with 


i;> 


l!c(i+1)  -C(1)||2  -  ||R<i+2>  -R(1+1)||2, 


we  can  conclude 


(3.8)  (R<1+2)-R(1+l))+(C(1+l>-c(1))-  0  as  i  — *  ‘ 

Subtracting  (3-7)  from  (3-8)  yields 

(R(i+2)  -R(i  +  1))  -  (R(1  +  1)  -R(i))  -*o  as  i— *». 


Thus,  for  a  sufficiently  large  Nq  and  fixed  n^,  we  can  keep 


(r 


(N0+l+i)  (Nq+1  ) 

10,'^0  10,^o 


)  ,  i  -  1 , 2  ,  •  •  •  ,n. 


arbitrarily  close  to 


(r(V1>_p‘V  , 

0’J0  o,Jo 

This,  however,  contradicts  (3.4)  and  (3-5)  both  being  true. 

Since  {R^n^}  and  lC^n^3  are  bounded,  there  must  exist 

(ni )  (nj) 

convergent  subsequences.  Suppose  R  1  — *  R  and  C  — »  C. 

Then,  in  light  of  (3-3), 


JnJ  (n,  )  (n,  ) 

G  =G+R  1  +C  1 


and 
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A(ni+1) 

G  X 


(n.+l)  (n  ) 
G  +  R  1  +  C  1 


both  converge  to 

£ 

G  = G+R+C  (in  anticipation  of  this  being  the  desired  solution) 

Since  8  is  an  element  of  and  is  an  element  of 

n »  ws  know  that  G*  £  fl  (these  cones  are 
closed).  Furthermore, 

( G-G*  , G* )  =  (G+R-G* ,G*)  -  (R,G*) 


=  lim(G+R 
i  -#oo 


(ni)  JV  ~(ni}  A^.+l)  i\Cnj  *^1)  \ 

-G  1  ,G  1  )  +  lim(G+C  1  -8  1  ,8  1  } 

i  -*co 


=  0+0. 


Similarly,  if  V  €  fl  Kc , 

(G-G* ,V)  =  (G+R-G*, V)  -  (R,V) 

(n  )  Jn  )  (n  )  (n  +1) 

=  lim( G+R  -G  ,V)  +  lim(G+C  1  -G  1  ,V)  s  0+0. 

i-»®  i  -*® 

Thus  G*  is  the  desired  solution  by  (2.1)  and  (2.2).  Moreover, 
since 

-C  minimizes  ||G+R-Fi|2,  F  €  K*w 

c 

and 


-R  minimizes  j|G+C-Fj|  ,  F  €  K*  , 


we  may  use  the  distance  reducing  property  of  projections  to  say 


||C(n)  -CH2  =  1|G  +C('n)-(G+C)|l2  St  iiR(n+1)-Ri|2 


=  ||G+R(n+1)  -  (G+R)  j|2  s  |ic<'n+1^  -  C  j| 2  for  all  n, 


Thus 


R^n^  — »R  and  — >C  as  n 


which  implies  that 


=  G+R^  +C^n_1^  and  =  G+R^  +  C^ 


both  converge  to 


G  =  G  +  R  +  C  as  n 


4.  OTHER  POINTS.  It  is  important  to  note  that  the  solution 
G*  =  G  + R  + C  does  not  uniquely  determine  R  and  C.  In  fact, 
if  we  begin  with  a  column  smoothing  rather  than  row  smoothing 
we  will  obtain  different  limiting  values  from  R  and  C  even 
though  the  same  limiting  G*  is  obtained. 


As  one  would  expect,  this  procedure  works  equally  well  when 
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the  order  restrictions  are  modified  to  require  nonincreasing 
rows,  or  nonincreasing  columns,  or  both.  One  has  only  to 
change  the  one  dimensional  smoothing  to  operate  in  the  appro¬ 
priate  direction.  The  procedure  can  also  be  adapted  to  higher 
dimensions,  although  in  this  case  the  number  of  required 
smoothings  quickly  becomes  large. 

We  also  wish  to  point  out  that  G*  itself  solves  many  more 
minimization  (maximization)  problems  than  the  least  squares 
problem  stated  above.  For  example,  from  Theorem  1.10  of  Barlow 
et  al.  ,  if  4  is  an  appropriate  convex  function  and  cp  is  a 
subgradient  (basically  a  derivative)  of  4,  then  G  solves 
the  problem 

a  b 

(4.1)  Maximize  E  E  { 4( f,  , )  +  (g,  ,-f,  , ) <p( f,  , )  ]w, , . 

F€K  n K  i=l  j  =  l  1J  J 

r  c  d 

Along  somewhat  similar  lines.  Theorem  3-1  of  Barlow  and 
Brunk  guarantees  that  the  problem 

a  b 

(4.2)  Minimize  E  E  ( * ( f, , )-g, , f, , >w. , 

F€K  D  K  i=l  J-l  J  J 

r  c  0 

—  1  # 

is  solved  by  (cp  where  once  again  i  is  an  appropriate 

convex  function  and  cp  is  a  subgradient  of  4. 

|| 

Thus  G  solves  a  much  wider  range  of  problems  than  is 
readily  apparent.  For  example,  suppose  one  has  a  multinomial 
random  vector  >  where  the  cell  probabilities  p^  are 

placed  in  a  rectangular  grid  and  one  wishes  to  find  the  maximum 
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likelihood  estimators  for  the  p^j  subject  to  nondecreasing 
(nonincreasing)  rows  and  columns.  This  problem  can  be  phrased 
in  terms  of  (4.1)  from  which  it  follows  that  the  solution  is 
given  by  G  where  G  =  (X^j/n)  and  =  1. 

Similarly,  if  the  are  independent  binomial  (n^.p^) 

random  variables,  one  can  show  that  the  maximum  likelihood  esti¬ 
mators  for  the  p^j  subject  to  nondecreasing  (nonincreasing) 
rows  and  columns  is  given  by  G*  where  G  =  (X^/n^j)  and  w^  =nij 

Finally,  in  order  to  illustrate  the  algorithm  on  a  larger 
table, we  considered  the  data  presented  in  Table  4.  The  entries 
are  the  first  year  grade  point  averages  of  2397  students  who 
entered  the  University  of  Iowa  in  the  Fall  of  1978.  The  inde¬ 
pendent  variables  are  the  composite  scores  on  the  ACT  Assess¬ 
ment  and  the  student's  high  school  percentile  rank.  The 
expected  first  year  grade  point  average  is  assumed  to  be  a  non¬ 
decreasing  function  of  both  of  these  Independent  variables. 

(The  number  in  parentheses  is  the  number  of  students  in  the 
category . ) 

The  least  squares  solution,  correct  to  four  significant 
digits,  was  obtained  after  500  Iterations  (250  row  smoothings 
and  250  column  smoothings)  at  a  cost  of  9  seconds  CPU  time  and 
84  cents.  These  results  are  given  in  Table  5  with  the  level 
sets  indicated.  Since  the  cost  of  our  algorithm  is  essentially 
linear  in  the  number  of  points  in  the  grid,  even  very  large 
arrays  can  be  isotonized  at  a  reasonable  cost. 
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